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Abstract 



Using the recently introduced Stripper approach to double-real radiation, we evaluate the 
total cross sections for the main partonic channels of the next-to-next-to leading order contri- 
butions to top quark pair production in hadronic collisions: gg — )■ ttgg^ gg — ^ ttqq^ qq — > tigg^ 
qq — >■ ttq'q' , q' ^ q. The results arc given as Laurent expansions in e, the parameter of dimensional 
' regularization, at a number of m/EcM values spreading the entire variation range, with m the 

mass of the top quark and Eqm the center-of-mass energy. We describe the details of our imple- 
i' mentation and demonstrate its main properties: pointwise convergence and efficiency. We also 

1^ . prove the cancellation of leading divergences after inclusion of the double- virtual and real- virtual 

contributions. On a more technical note, we extended the double-soft current formulae to the case 
of massive partons. 



1. Introduction 



> 
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■ In we have proposed Stripper (SecToR Improved Phase sPacE for real Radiation), a 
novel subtraction scheme for the evaluation of double-real radiation contributions in next-to-next- 
to leading order (NNLO) QCD calculations. As suggested by the name, it is the phase space 
that acquires a special role in our approach. Once it is suitably parameterized and decomposed, 
Laurent expansions of arbitrary infrared safe observables can be obtained without any analytic 
integration, by applying numerical Monte Carlo methods. An important feature of Stripper is 

... its process independence. In the actual calculation, general subtraction terms are combined with 

r> I process dependent amplitudes. The simplicity of our construction contrasted with the complexity 

■ of double-real radiation singularity structure may lead to scepticism. The present publication is 
meant to prove that Stripper delivers on its promises. 

There are many other subtraction schemes for real radiation. At next-to- leading order (NLO), 
most calculations are done with the method of 0, 01 , but other approaches are subject to active 
development 0, 01 (see also @). At the NNLO level, the situation is more complex. Much has 
been achieved with Sector Decomposition [7|, |8|, |9| and Antenna Subtraction 

Einilillil, but 

there are more specialised methods for colorless state s 1141 and, very recently, for massive final 



states [15[. General tools are also being developed in [16|, [17|, [18|, [19|, |20[. There were, of course, 
many other proposals 21, 22, 2^, which have not been completely developed. 

To demonstrate the virtues of Stripper, we have chosen hadronic top quark pair production, 
as it is of great phenomenological relevance, but its theoretical description is still incomplete. The 
current state of the art in the field is as follows 

1. differential cross sections including complete off-shell effects and leptonic decays are known 
at NLO [3, [25! (until recently, they were only available in the narrow- width approximation 
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fixed order threshold expansions for the quark;-annihilation and gluon-fusion channels are 
known at NNLO up to constants, both for the total cross section ^2^], and for the invariant 
mass distribution 29]; 

soft-gluon resummation for the previously mentioned channels is understood at the NNLL 



level for the total cross section [30|,|3l[ (see also |32|), and selected differential distributions 
^ (see also jsij); 

4. mixed soft-gluon and Coulomb resummation is understood at the NNLL level as well [ssj ]: 



5. two- loop virtual amplitudes are known analytically in the high-energy limit 36l. |37|. for the 
planar contributions 3^ 3^, and fermionic contributions in the quark-annihilation channel 
[40|; in the same channel, the full amplitude is known numerically 41 1; 

6. one-loop squared amplitudes are known analytically at the NNLO level [H, |4§|; 

7. one-loop real-virtual (one additional massless parton in the final state) amplitudes have 
been obtained in the course of several projects connected to top quark pair production in 
association with jets 0, |4^ |4^ ; 

8. approximations of the one- loop real- virtual amplitudes are known in all singular limits in- 



volving only massless partons [47|, |48|, |4^ i51|| (this is needed for the evaluation of the 
real-virtual contributions). 

We are interested in an NNLO calculation. Schematically, the involved partonic cross sections 
are a sum of three terms mentioned already above 

daf^f = da""^ + da'''' + da'''' , (1) 

where da^^ denotes the double-virtual (two- loop and one- loop squared), da"^ the real- virtual 
(one- loop with one additional parton), and da"" the double-real (tree-level with two additional 
partons) corrections. We will ignore the need for collinear renormalization, which involves lower 
order cross sections expanded to higher orders in e, the parameter of dimensional regularization. 
These have been derived (although not published) for the analysis of [H^ . 

Currently missing are the double-real and real-virtual contributions. Once these are known, 
one can provide complete NNLO cross sections beyond the known threshold expansions. The 
real-virtual contribution has a simpler singularity structure than the double-real, and should be 
obtainable once the soft-gluon current in the presence of heavy quarks has been derived. 

As far as the double-real contribution is concerned, there are many partonic channels, which 
need to be considered in principle. Nevertheless, current phenomenological applications require 
the knowledge of the cross sections with gluon-gluon and quark-anti-quark initial states. This 
leads to our choice of the following channels 

99 ttgg, gg ttqq, qq ttgg, qq ttq'q' . (2) 

We will only consider the case q' ^ q, because the case of identical quarks is expected to be 
numerically irrelevant. In the present publication, we will provide numerical values for the cross 
sections for the four channels as function of m/EcM- 

The paper is organized as follows. In the next Section, we will discuss the phase space, its 
volume, parameterization and decomposition. Subsequently, we will describe the derivation of 
the subtraction terms, convergence and cancellation of leading divergences. In Section 4, we will 
describe the technical details of our implementation, demonstrate its efficiency and describe several 
tests. Section 5 contains the results of numerical simulations for the four chosen channels. Apart 
from the main text and Conclusions, the publication consists of a number of Appendices. They 
contain the collinear splitting functions, a discussion of the double-soft limit in the presence of 
massive partons, another approach to the collinear limit in the double-soft limit, the Born cross 
sections, and finally a list of the software that we have used. 
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2. Phase space 

2.1. Volume 

A numerical approach, as the one advocated here, poses substantial problems, when assessing 
the correctness of both the approach and implementation. We will, therefore, start by introducing 
the only truly non-trivial integral that relates to our computation, but can be evaluated entirely 
analytically: the volume of the phase space. 

We are interested in the following class of processes 

a{pi) + b{p2) t{qi) + i(q2) + c(fci) + ^(fcs) , (3) 

where a, b are initial and c, d final state partons, in particular ah — gg or qq, and cd — gg or q' q' . 
Additionally 

pl^pl = kl^kl = Q , ql = ql = m V , (4) 

and as usual 

S={pi+P2f. (5) 

Throughout this publication we will work in the partonic center-of-mass system. 
The phase space measure in c?-dimensions is 

d'^-^ki d'^-^k2 d'^-^qi d'^-^q2 ^rfe■WW, , n 
= {2.Y-^2kl i2nr-^2kO (2.)'^-i2<z? (27r)''-i2gO ('"^ Hk, + k2 + q^ + q2 - P, - P2) , (6) 

with d = 4 — 2e as usual. We wish to evaluate the integral of unity with this measure. The result 
will be given as a product of two functions 

J d<l>i = Pi{s,eMx,e) , (7) 

with 



1-/3 ^ / 4m2 

and P4 the volume of the phase space in the purely massless case 

P4(s,e) = 2^"+6^^-5+3^ ^ rs'-'^ , (9) 

^ ' ' r(3-3e)r(4-4e) ' ^ ^ 

which can be easily obtained from the imaginary part of the three-loop massless sunrise diagram 
[sij . By definition, the $ function must satisfy two boundary conditions 

$(l,e) = 0, $(0,e) = l. (10) 

The first of the above equations is just the vanishing of the phase space at threshold, which is 
located at x = 1. The second follows from the normalization in the massless case, which in turn 
corresponds to a; = 0. 

In order to obtain <I>, we will use the method of differential equations 5J, |55|- To this end, we 
exploit again the fact that $ is given, up to normalization, by the imaginary part of the three-loop 
sunrise diagram, this time with two massless and two massive lines, see Fig.[TJ Reducing the mass 
derivatives of the three occurring master integrals using integration-by-parts identities, we obtain 
the following two equations 

^$(a;,e) =^ ^ 2(x - 1) e e) + (1 -|- x)^ o h (H) 



dx \ x^ ) x^ \ dx 

d^-^ix.e) l-2Qx + ?,x^ d-^ix.e) 24 ^, ^ 
i . i i . L "^{x e) 

dx'^ a;(l — a;^) dx x{l -\- x)"^ ' 

'l-22x + x^ d^{x,e) 24(2 - e) 2(1 - 2e 



x{l — x"^) dx x{\ + xY {1 + x) 



M'(a;,e) + -^— /$(a;,e) , (12 



3 



Pi +P2 




I 



Figure 1: Three- loop sunrise graph used to obtain the phase space volume. Thick lines are massive, whereas the 
dashed line denotes a cut, which in this case corresponds to the imaginary part of the integral. 

where we have only kept two of the master integrals, $ and Notice that ^f's sole purpose is to 
provide a solvable second order differential equation and a neat boundary condition, which follows 
again from the vanishing of the phase space at threshold 

*(l,e) = 0. (13) 

For this reason, we do not even bother specifying its exact definition. 

With the three boundary conditions, the solution of the system of differential equations is 
unique, and can be obtained recursively as a series expansion in e. In principle, we need five 
terms of the expansion corresponding to the five relevant terms of the expansion of the cross 
section, ranging from l/e"* to e°. As the expressions quickly become extremely lengthy, we will 
only reproduce the first two, and give a high precision numerical value at some benchmark point 
for the complete expansior0. $ is given by 

oo 

6'$W(^) ^ (14) 

i=0 



The analytic result for $(a;) is attached to the electronic preprint version of this publication 
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with 



48a;2iJ(-l,0,a:) , 24:X^H{0,0,x) , 12 {x^ + 5x^ + 6x^ + 5x + l) xH{0,x) 

(15) 



(a; + 1)4 ^ (a; + 1)4 ^ (x + 1)6 

f 23x4 + (34 + 47r2) + [iir^ - 34) a;^ - 23a; - 1 



(a: + 1)5 

672a;2i?(-l,-l,0,a;) 96x2iJ(-l, 0, -1, x) 144x2i/(-l, 0, 0, x 



(X + 1)4 (X + 1)4 (X + 1)4 

480x2ff(-l,0, l,x) 336x^^(0, -l,0,x) 48x2i7(0, 0, -1, x) 

(X + 1)4 + (X + 1)4 + (X + 1)4 

72x2i7(0,0,0,x) 240x2i/(0,0,l,a;) 4 (Sx^ + 56x + 3) xi/(0, 0, 



+ 



(X + 1)4 (X + 1)4 (X+1)4 

8 (15x4 + 122x3 + 190x2 + 122x + 15) xH{-l, 0, x) 

(X + 1)6 

24 (x'' + 5x3 ^ g2;2 + 5a; + 1) xH{0, -1, x) 

(x + l)6 

120 (x4 + 5x3 _,_ g^2 + 5a; + 1) xH{0, 1, x) 
(x + l)6 

4 (2x5 + 61x4 + 4 (50 + 7r2) x^ + 2 (89 + 47r2) x^ + (86 + 47r2) x + 13) xi7(0, x) 

(x + l)6 

(-2x5 - 46x* + 4 (87r2 - 17) x^ + 4 (l7 + 87r2) x^ + 46x + 2) H{-1, x) 

(X + 1)5 

10 (x5 + 23x4 + 34x3 - 34x2 _ 23x - l) H{1, x) 

(x + 1)5 

2x (tt^ (18x4 + 43x3 + Sx^ + 43x + 18) + 180x(x + 1)^(3) 



3(x + l)e 



(16) 



The H functions are standard harmonic polylogarithms (HPL) [56j. Their initial weight is 2, but 
the last term of the expansion we are interested in, i.e. e'^, contains HPLs up to weight six. 

We can also obtain the behavior of $ near threshold, either from the differential equations, or 
from the actual solution. The result is 



i:p|,^)=^/3^-^°^(l + 0(e))+0(n- (17) 



The numerical benchmark expansion is chosen at x = 1/2 

4>(l/2,e) = 0.00001122829901964763 

+ 0.0001283543727784153 e 

+ 0.0007325963156455679 

+ 0.002782712436211506 £3 

+ 0.007910064621069109 e"* 

+ 0{e^) . (18) 

One may wonder why it is not sufficient to test the implementation close to the massless case 
and avoid the work that led to the above expressions. The reason is that the presence of large 
logarithms of the mass, up to log"* (m^/s), implies a high sensitivity of the result to the value of the 
small mass. A point like x = 1/2 has no special properties and, therefore, no large cancellations 
are expected. 
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Using the result for and its two derivatives in x, we have a complete set of master integrals 
and can evaluate the integral of any polynomial in scalar products of pi +P2, ^1,2 and qi^2- We 
will find it later useful to have the result for the following integral 



'^'^^ \ s J " 48(e-l)2(12e2-31e + 20)(l-x)(l + a;)4 

X ((a; + 1) {x^ - 1) x^ (e^ {6x^ - 8a; + 6) + e {-19x^ + 4a; - 19) + 2 {7x^ + 3x + 7)) $"(a;) 

-{x + l)x{2e^ {5x^ - 66x^ + S2x'^ - 66a; + 5) + {-47x^ + 526x^ - 2\^x^ + 510a; - 35) 
+e (72a;'' - 636a;3 - 66a;2 - 628a; + 34) + 4 (-ga;"* + 59a;3 + 26x2 + 62a; - 2) ) $'(a;) 
-2(a; - 1) (2e'* (a;'' + 60a;3 + 26a;2 + 60x + l) - (lla;'' + 562a;3 + 470a;2 + 562x + ll) 
+6^ (22a;'* + 933a;3 + 10283;^ + 933a; + 22) - e {\9x^ + 660a;3 + 8523;^ + 660a; + 19) 

+2 (3x'* + 85a;3 + 122^^ + 85a; + 3) ) $(x)) . (19) 



The value at our benchmark point is 

(1.4553533 + 16.673868 e + 95.377076 <? + 363.03219 + 1033.8027 + 0(e^)) x 10"^ . 

2.2. Parameterization oj the massless system 

We will now introduce a suitable parameterization of the phase space. We will closely follow 
the lines of [ij, where the massless system has been specified. In the next subsection, we will 
define a parameterization for the heavy system. 

Before we give the momentum vectors, let us note that we can always choose them such that 
their e-dimensional components vanish. This is due to the rotational invariance remaining in the 
system as long as we only have three vectors Pi,fci,fc2 (notice that p2 = —pi by assumption). 
Therefore, we will specify the vectors, as if they were purely four-dimensional. The only conse- 
quence of the existence of the additional degrees of freedom is the modified form of the phase 
space measure, which is then sufficient to regulate all singularities. We will also exploit rotational 
invariance and space inversion invariance of the matrix elements (which can also be viewed as 
d-dimensional rotation invariance) to restrict the momenta as follows 

fc^ = , k^>0. (21) 

In the actual Monte Carlo simulation, the momenta should be rotated randomly around the z-axis 
and the sign of a;-axis should also be chosen at random, in order to fill out the complete phase 
space. 

With the above assumptions, let 

p'l = ^^(1, 0,0,1), 

P^2 = ^(1,0,0,-1), 

< = ^/3'(l,O,sin0i,cos0i) , 

/32(l,sin</)sin6'2,cos(/)sin^2,cos02), 



2 

ki = ll , 

k!^ = (22) 
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where 0,^1.2 G [0,7r], and auxiliary vectors needed to define soft limits, whereas are 

used to parameterize the energies. Notice the hats above the variables. We shall use them to 
denote all variables, which are going to be transformed due to further phase space decomposition. 
The angular variables are replaced by another set in two steps. 
We first define 

'71.2 = ^(1 -cos 01,2) , (23) 

V3 = ^(1- COS 6*3) 

= -(1 — cos0sin(?i sin6'2 ~ COS01 cos(?2) 

= i(l-cos(6'i-6l2) + (l-cos0)sin6lisin6'2) , (24) 

where ^3 is the relative angle between ki and fc2, and by definition r)i, 2,^73 G [0, 1]. One of the 
main ideas of the subtraction scheme [H is to change variables in such a way that all collinear 
limits be parameterized with just two variables, fji and 772. In order to do so, we introduce 

. 1 (1-COS(01-02))(1 + COS0) 

^ 2 1-cos(6li-6l2) + (l-cos</))sin6lisin6l2 ^ ^ 

which can be inverted to give 

m ■ (26) 

7)1 + 7)2 - 2fiifj2 - 2(1 - 2C)v??i(i - m)m{^ - m) 

Clearly, the collinear limits are now at 771 = 0, 7)2 = 0, r}i = 1, 7)2 = 1 or 771 = 772. While 61^2 are 
obtained from Eq. (^5]) . (j) is given by solving Eq. (IM)) and Eq. (1^ 



l-(l-277i)(l-2^2)- 



, »)l+»)2-2))ir)2-2(l-2C)^/(l-i)i)ih(l-»)2)»)2 
C0S(/)= = . (27) 

4v(i - '?i)??i(i - m)fi2 

Notice that 

??i = ?72 cos = 1 , (28) 

whereas 

7)1 = V 772 = V 771 = 1 V ?72 = 1 ^ cos = 2C - 1 . (29) 

The last statement is valid, when 771 ^772, but seems to contradict implication (j28l) . Fortunately, 
in the limiting cases ^1 = 772 = and 771 = 772 = 1 the momentum vectors do not depend on 0. 
The final set of parameters specifying the kinematics of the massless partons is 

C, m, fl2, ii, 6 • (30) 

The first three are unrestricted within the range [0, 1], whereas the energy variables belong to one 
of the two non-overlapping regions (apart from a measure zero set) [l| 

{ (li, 6) : < a < 1 , < I2 < a UaAii)} , (31) 

{ (li, 6) : < I2 < 1 , < li < 6 ?™ax(6)} , (32) 

where 

a_(0=^min(^l, i^^^)<l. (33) 

Not only do these conditions guarantee that the massive states can always be produced, but they 
are also suggestive of a decomposition of the phase space, which we will perform later on. 
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Having specified the parameterization of the phase space, we can rewrite the measure Eq. ([B]) 
in the new variables. We split it into two parts 



(i$4 = d*3(Pi +P2\ki,k2) d$2(Q; '71,92) , (34) 

with 

Q = Pi + P2 - fci - ^2 • (35) 

d^s,{pi +P2; fci, ^2) is not exactly the three-particle phase space of fci, k2 and Q, because the only 
constraint that it is subjected to is > Arn^. On the other hand, d$2(Q; 91, 92) is the two-particle 
phase space. We have 

d<^,{p,^P2-MM = 8(2^)5r(i - 26) "''''^'"'' ^^^^ ~ 

X if,^{l - 77i))-(ry2(l - ,}2))--^4^ 

X dC dfji dfj2 d^i d^2 ■ (36) 

The first line above will be of no further concern, since we are only going to perform variable 
changes on the subset 571,572,^1,^2- Therefore, we will define 

dy^^i = (r?i(i - 7?i))--(7?2(i - m))-' il'^'il'^' dm dm dii di^ , (38) 

with c?4>3 = d^j.(; dfj,,i^. Despite the splitting, d/i,,^ depends on ( through 773. 

2.3. Parameterization of the massive system 

In order to parameterize the massive system, we perform a boost to the center-of-mass frame 
of Q. Denoting the momenta of the heavy quarks in this frame by q[ 21 '^^ have 

9° = 




qt ^ ql + I q-> + —^-11= ] , ^ ^ 1,2 . (39) 



The problem that we face is that once three {d~ l)-dimensional momenta of the massless partons 
have been specified with e-dimensional components vanishing, we do not have the freedom to 
keep the latter components of q'l 2 vanishing anymore. The easiest solution would be to restrict 
the momenta of the heavy quarks, which are always resolved after all, to lie in the four physical 
dimensions. This would simplify the parameterization, but we would loose the possibility to use 
the integrals derived in Section 12.11 Furthermore, to obtain finite partonic cross sections, it is 
necessary to add collinear counterterms, which are convolutions of splitting functions with lower 
order cross sections. If we would like to use the results obtained in [52] for this purpose, we need 
the heavy quarks in d-dimensions. 

Let us, therefore, define the q[ 2 momenta through three spherical angles 9q, (pq and pq. 



<7i = -92 = \VQ^TW~l (40) 



— -92 — 2V^"^/-'" — 

(sin pq sin (j)Q sin 9q n^'^"'*-' , cos pq sin 0q sin 9q , cos (j)Q sin 6q , cos 9q ) 
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In principle, the three angles should lie in the range [0, tt]. Nevertheless, we can assume (j)Q G [0, 27r] 
and pq e [0,7r/2], as long as we exploit the independence of the results from the sign of n^'^~^\ 
In fact, without loss of generality, we can set 



(41) 



and forget about the (d — 5)-dimensional components. Thus we have to work with five-dimensional 
vectors. We will soon see that the contribution of those vectors, which have a non-vanishing fifth 
dimension is suppressed by a power of e as one would expect. 
The two-particle phase space is now 



rf$2(Q;'?i,g2) 



8(27r)2r(l - 2e) 
4i+T(-2e) 

r2(-e)(l-cos2pQ)'+^ 
fi2 d cos 9q d(j)Q d cos pq 



Am? 



l-2e 



(1-COS^ 



sm 



d cos 9q d(f)Q d cos Pq 



(42) 



It depends on C, 771,2, ^1,2 only through , although the momentum vectors gi,2 depend on each of 
these variables independently. Close to threshold, where sa s, we recover the behavior Eq. (|17p 



2-3e o9-10e 



d$4 OC S^'-'^P 



(43) 



More interestingly, however, the ratio r{—2e)/r'^{—e) is of the order e, which means that we need 
a divergent contribution from the integral to obtain a cross section in four dimensions. This is 
indeed guaranteed by the following 



4i+T(-2e) 



r2(-e)(l-cos2pQ)^ 
where the "4-" -distribution is defined as 
1 



6{l - cos pq) + 



4i+^r(-2e) 



1 



(1 - C0S2 Pq) 



1+e 



(44) 



dcos 



PQ 



(1-C0S2 Pq)1 + ^ 



/(cOSpg)-^ dcOSpQ-^^— ^^1— y^(/(cOSpQ)-/(l)) , 



(45) 

and the integrand on the right-hand side should be expanded in a Taylor series in e. While we leave 
the discussion of the implementation details to Section |4l we note that we chose to use equation 
Eq. (j44|) to divide the phase space into two contributions 



P2 "Cost 



p!"2^dcos 9q 



d cos Pq 



(46) 



with 



M2 



(4^)T(1 - e) 
8(27r)2r(l - 2e) 

4(27r)2r(-e) 
1 



(1-C0S2 Pq)1+^ 



'1- 



4to2 



(1 - cos^ 



^-^(sin2</>Q)- 



(47) 



'1- 



4to2 



(1 - cos^ 



(sin^ <PqY 



(48) 



d^''2^'^^ would be the entire phase space, if we could rotate the e-dimensional components away. 
One can expect that the additional contribution from d'^^2' '^iH be small in practice. We will show 
later that this is indeed the case. 
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At this point, we would like to note that the adopted solution to the problem of a c?-dimensional 
phase space for the heavy quarks is by no means unique. One could, for example, use the fact 
that the e-dimensional components of the heavy quark momentum vectors are only relevant to 
the terms singular in e, which are obtained after one of the massless vectors has been removed (at 
least one soft or coUinear limit). We then have only two {d — l)-dimensional vectors, and could 
rotate away the spurious components of qi^2- This approach would only be correct, if the reference 
frame for the parameterization of gi^2 were defined in relation to pi and fci + fc2. This in turn, 
would be a simplification for the massive system, but a complication to the decomposition of the 
phase space, which we want to perform in Section 12.41 

2.4- Decomposition 




(S9 

Figure 2: Decomposition of the phase space in the triple-collinear sector. The variable substitutions, which map 
the integration range onto the unit hypercube are specified. Furthermore, ^2 = ^maxi^i) and the second branch 
starting with the dashed line is symmetric to the first. 

The last step of our treatment of the phase space is a two-level decomposition according to 
singularities. At the first level, we partition the phase space with suitable selector functions. The 
latter are defined on the phase space, add up to unity, and regulate part of the divergences. In 
particular, we introduce a selector function for the triple-collinear sector, in which we allow for 
coUinear divergences due to partons with momenta pi, fci and /c2, but not p2- There is also a 
symmetric function that does just the same upon replacement of pi with p2, but we ignore it, as 
its contribution can be recovered without additional computation (see Section |4l). Moreover, we 
introduce a selector, which allows for coUinear divergences due to ki being parallel to pi, or fc2 
parallel to p2, but no other configuration. This function defines the double-collinear sector, and 
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6 > f 1 



Figure 3: Decomposition of the phase space in the double-coUinear sector. The notation is as in Fig. [21 



has a symmetric counterpart in pi P2, which we again do not discuss any further. The triple- 
and double-colhnear sectors may be overlapping in the sense that several selector functions do 
not vanish for a given momentum configuration. The only condition is that the divergences are 
properly regulated. In we have given two examples of selector functions, which achieve this 
goal, one for our present problem, and one completely general for any number of massless final 
states. Apart from numerical efficiency, nothing depends on the choice. In the present work, we 
define the sectors implicitly as follows 

1. n\> {) A n| > —a ^2 in the triple-collinear sector; 

2. nf > A < "2 the double-collinear sector, 

where a > is an arbitrary parameter, which we will take to be a = 1/2 (we checked independence 
of some results on this parameter). Notice that the sharp cut > in both cases is necessary for 
the later use of symmetries. Moreover, we have defined the conditions through the vectors 
because in the strict soft limits the actual momentum vectors vanish and it is impossible to check 
to which sector they belong. 

Having simplified the problem as far as the type of coUinear singularities is concerned, we 
will perform a second level decomposition. The purpose is to factorize the divergences of the 
propagators in the amplitudes. The set of offending invariants is 



Sl5 


- bi 




= -sp'^iifji , 




•SI6 


= {Pi 


-k^f 






S26 




-k,f 






S56 


= ih 


+ k2f 


= sf3^iii2ri3 , 




■Sl56 


- {Pi 




k2f = -sP^iim - 


^-(2^2 - f3^iii2V3) , 


S256 


= {P2 




fc2)' =-5/3^(6(1 - 


' Til) + ' m) - P''iii2m) 



Assuming that we are only concerned by the coUinear singularities, sis, sig, S56, S156 are relevant 
to the triple-collinear sector, whereas si5,S26 to the double-collinear sector. In case partons may 
also become soft, there will be soft-coUinear singularities in the double-collinear sector due to S156 
and S256- Purely soft singularities may involve other propagators, in particular the massive, with 
a general form as follows 

{p + hf-p^ = 2i,p-n,, i^l,2, 
{p + ki + k2f-p^ = p ■ ni + £,2 P ■ n2 + iii2 ni ■ n2) . (50) 
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Due to the polynomial form of the singular denominators, it is possible to factorize the divergences 
with a change of variables in such a way, that each expression be a product of variables to some 
powers, the variables themselves inducing divergences when vanishing, and a regular function. 
This is the well known sector decomposition Q. The difference to the usual treatment is that, due 
to the process independent nature of singularities in QCD, we do not have to consider abstract 
expressions, but can view the decomposition as a choice of the order of the soft and coUinear limits. 
A complete schematic representation of the variable transformations, which lead to factorization 
of divergences and specify the order of singular limits is given in Fig. [2] for the triple- and Fig. [3] 
for the double-collinear sector. Each level in these decomposition trees factorizes a certain type of 
singularities. For Fig. [2l we have 

I) factorization of the soft singularities; 

II, III) factorization of the collinear singularities; 

IV) factorization of the soft-collinear singularities. 

The case of Fig. [3] is even simpler, as some of the levels disappear. Note that for each tree, we 
start with a new set of variables 7yi_2 = Tji,2 and S.1.2 = £.1.2 at the root, and continue with the 
substitutions down to the leaves. The missing right branches corresponding to ^2 > £.1 can be 
recovered by changing the order of the final state partons. For this reason, we will ignore them 
as well. Finally, the variable substitutions guarantee that 771,2,^1,2 € [0, 1], as the range Eq. (|3ip 
gets expanded. 

Proving that this procedure is sufficient for factorization is a simple matter of substitutions 
and can be done with pen and paper. We shall not reproduce the uninteresting transformations 
here, since we shall not need them anymore. Nevertheless, we point out that the factorization 
of ?73, present in 555, introduces two regular functions, which we will encounter later. They are 
defined as 

f?3i(??i,?72) — V3 (m,7;mV2] (51) 



(2-^2)^ 



(2 + 7^2(1 - 2771) - 2(1 - 2C)v/r?2(l-77l)(2- 7/1772)) ' 



1 

,. ...J , ..J., —'I 

mvi V 2 



7732(771,772) = 773 ( 771, -7/1(2 - 772) ) (52) 



2 (2 + (1 - 27/i)(2 - 7/2) - 2(1 - 2C)v/(l-m)(2 - 772)(2 - 771(2 -772)) 

with 773(7)1,7)2) as in Eq. (|26|) . The first of these functions is relevant to the sectors S(,...,S^, 
whereas the second to Sl,S[. The derivation of the subtraction terms in Section IXTl requires the 
knowledge of the behavior of 7/31(771,772) for small 7/2 



r?3i(?7i,??2) = 1 + (1-20^2(1 -m)'72+ 0(7/2) . (53) 

Tab. [1] contains the complete set of variable transformations, i.e. the expressions of the original 
kinematic variables 7/i,2:Ci,2 in terms of the sector variables 7712,^1,2- For the double-collinear 
sector, we have used the transformation 772 ^ 1 — 772 first, in order to profit from the phase space 
measure Eq. (j38p . which is invariant under this transformation. 

The above factorization procedure does not leave the phase space measure unchanged, even 
though it is only the factor dii^i^ that is transformed. The latter assumes, for a given sector S, 
the following form 

dM,? = ^a.+b.e^a2+b,e^a,+b,.^a,+b,. ^^eg drJ,d7^2d^^d^2 , (54) 

where the first factor regulates the divergences, whereas the second is a regular function, which 
we will later on expand in e in a straightforward Taylor series. Both factors are given in Tab. [2j 
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Table 1: Original kinematic variables, 171 , »?2i €l i ?2! expressed through the sector variables, »?i , »72i €i > 52i with 

§2 = (max{S,l)- 

2. 5. Normalization of the cross section 

Although we have now specified the four-particle phase space completely, there is still a freedom 
in defining the divergent cross section for double-real radiation by including a function, which 
equals unity, when e = 0. Usually such a function is introduced due to the e-dependence of the 
bare coupling constant, which, in the MS scheme, is chosen to be 

«° = (^) Z^Mf^'),^) , (55) 

where fi is the renormalization scale and Za^ the renormalization constant. Of course, in our 
case Za^ — 1, since we are integrating tree-level amplitudes. The problem with this approach 
is that our matrix elements are proportional to a^, whereas the dimension of the phase space is 
the same as that of a three-loop integral. This would leave unbalanced, large and dimensionful 
factors of /i^*^. We can compensate them by multiplying the cross section with the inverse of the 
parenthesized factor in Eq. (j55p to power e. Equivalently, we will express all tree level amplitudes 
by the renormalized strong coupling and use the following definition of the partonic cross section 

J d$4 Fj , (56) 

where the overline over the matrix element squared signifies the sums and averages in color and 
spin, as well as statistical factors for identical final states. Fj is a jet function defining the 
observable O. In what follows, we will mostly use a trivial jet function Fj = 1. Notice, however, 
that a test in Section |4] will be performed with a non-trivial Fj. 

Let us stress, that we could have considered a different factor in Eq. (fSG)) . as long as we would 
include it in all the other contributions that enter the finite physical cross section. From the point 
of view of physics, this factor is irrelevant, and yet it can play a substantial role in obtaining precise 
numerical values. For example, multiplying by Z?^"*^, we would remove dominant logarithms of /3 
close to threshold, which would substantially lower the double-real cross section, while enhancing 
the others. We can also lower the contribution of the finite part of the cross section by multiplying 
with a constant, but e dependent factor. The decision on what to enhance and what to diminish 
can only be taken once all terms are implemented in a numerical program, because only then 
can we check, what contributes the largest absolute errors. We will leave this problem to future 
studies. 



1 //x e'^'' \ 
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regulator 


reg 
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Si 




((i-m)(2-m.2))-c7-^^(^^;^^^'^^)) 


1- 


-2£ 


Si 




((l-.2)(2-.1.2))-C7-^^('^^f^''^^)) 


1- 


-2£ 


Si 




((1 7?2)(2 mV2^2)) (2 ( 2-771^2 j 


l-2£ 


Si 


^l-2e„l-2c t3-4etl-2£ 
Vl V2 «1 ?2 


((1 - 77i)(2 - r?2)(2 - ,71(2 - V2))r^i-''vl2'%VuV2) 


Si 


l-2£ l-2£f3-4£(^l-2£ 
'/I '12 ?1 ?2 


((1 - r72)(2 - r7i)(2 - 772(2 - ViW ^i~''vl2'%V2,Vi) 


si' 


^r^^%^er'^e2^-^^ 


((l-r7i)(l-r72))-^^7-'^ 


( V3 ^ 

- '71 - m\J 


1^ 


~2£ 


si' 


vrv2'e^-''e2~'' 


i{l-V2)il-ViC2)r^2-'' 


( V3 ] 

V|i-'?2-??i6lj 


l-2£ 



Table 2: Integration measure, d^^j expressed through the sector variables, '»?i,2i^i,2, decomposed into the product 
of their powers used to regulate the divergences, and a regular function, f/g^, which can be expanded in e. 

3. Subtraction and integrated subtraction terms 

3.1. Derivation 

The decomposition of the phase space introduced in the previous sections is sufficient to derive 
Laurent expansions of arbitrary infrared safe observables. In order to obtain explicit expressions 
for a given sector S, we define 

Ms - vl^'''vl^'''ei^'''e2^''' JM^ , (57) 

where the constants have been defined in Eq. ([5^ . and are given for each sector in Tab. [2] The 
averaged matrix element, |A^4P, has been introduced in Eq. (|56|) . fffts must be regular, by infrared 
power counting in QCD, in limits of any of 771,2,^1.2 vanishing. This can be checked explicitly, 
with the formulae introduced later in this section. 
The cross section is now 

'^o=E'^o^ (58) 
s 
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where 



1 /^2g7E 

2s V 47r 



3€ 



rf^C c^^/i dm d^i d^2 d^2 Fj — 



d( drji dr]2 d^i d^2 dcost 



d cos pq h}rj , 



with the integrand 



.(5) 



1 

2^ 



47r 



3e 



MC Ms ^ 6*5 Fj 



l-6ie 
^1 



l-62£ fl-bse cl-b4e 
'12 ?1 ?2 



■ 



Ms 

(59) 
(60) 



The bi constants have been defined in Eq. (|54p . and are given for each sector in Tab. [21 The jet 
function Fj has been introduced in Eq. ((56|). Finally, 6*5 is the selector function described at the 
beginning of Section 12.41 We remind the reader that the full phase space is covered by changing 
the order of the final state massless partons, and swapping pi and p2, which can also be thought 
of as changing the order of the initial state massless partons. 

The Laurent expansion of the cross section contribution, ctq , is obtained by using 



1 



il-be 



b e 



00 

E 

n=0 



[bey 



ln"(A) 



where A = '71.2,^1,2, and the 



dX 



-distribution is 



ln"(A) 



/(A)= r^^(/(A)-/(0)) 
+ Jo 



A more practical, albeit equivalent, application of these formulae is 



dX 

Xi-b, 



/(A) 



dX 



/(O) /(A)-/(0) 
be Ai-''^ 



(61) 



(62) 



(63) 



In any case, Eq. ([5^)) involves four singular integrations and the above formula has to be applied 
iteratively. The result contains the integrand at sixteen different points obtained by setting the 
variables in all possible subsets of {r]i , 7]2 , £,1 , ^2} to zero. Using Eq. (155)) gives a convergent 
integrand 

(64) 



Considered differently, the terms in proportional to negative powers of e are called inte- 

grated subtraction terms, those free of singularities are simply called subtraction terms if any of 
the variables vanishes, just as in the classic approach to subtraction schemes. Notice that the 
only analytic integration that is needed here is the rather trivial integral of l/X^^'"^. This is the 
main difference to the traditional approach. The limits of c?<i>2, ^5^1 selector and jet functions 
are obtained by directly setting variables to zero, and are process independent. The only process 
dependent information is in QHg. The vanishing of the sector variables corresponds, however, 
to singular limits of QCD amplitudes. Thus, we can obtain the subtraction and integrated sub- 
traction terms from the splitting functions and soft currents exactly as it is done at NLO. The 
sector decomposition of Section 12.41 guarantees the independence of the result from the order, in 
which the limits are taken. The process dependent information will now be shifted to reduced 
d-dimensional matrix elements. 

In order to derive the relevant formulae, let X C {771,772,^1,^2} be the subset of vanishing 
variables in a given limit, and define 



lim Ms 

X-J-O 



g^{M3\Y\M3) or 



lim Ms 



/ {M2\Y\M2 



(65) 
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depending on whether the Born matrix elements above correspond to reduced processes with one, 
jA^a), or two partons, |A^2), less. V is an operator in spin and color space, and depends on 
the flavors of the partons involved. Let us also define a shorthand notation for the divergence 
regulating product of integration variables 



(66) 



We now consider the various limiting cases starting from those relevant to the triple-collinear 
sector. With the flavor assignment 



ai(pi) + 02(^2) t{qi) + t{q2) + 05(^1) + 06(^2) , 
there are nine cases, which are identified from X with the help of Tab. [T] 
1- '?i = j)2 = 



= lim ^s- 



(67) 



(68) 
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The splitting functions Pa^a^ae given in [Appendix A[ They depend on the following 
variables 

= x,=(3H,, xe=f3H2, (69) 

and 
with 



kl, = (0,0,1,0), 



(70) 
(71) 



m + m - 2{i - 2(:)Vrnm 

X (0, 2|77i- 772|VC(l-C), 2V^-(??i+'?2)(l-2C), 0) . (72) 

The last vector is symmetric and homogeneous in f)i and 772, fcj^g (^1,^72) = k^Q{fl2,'>)i) — 
KeC^^rii/m) = fc'|^6(l,??2/7?i). Moreover 



±5 



1 1-C 



(73) 



with 



fc^[o = (0,l,0,0). (74) 

This asymptotic behavior is necessary for Wis not to be singular. This limit is usually 
responsible for eight of the fifteen subtraction terms. 



m = 



2 P^^ 
lim ^g^lBL 



with 



1 - P^^i 



k^^ - (0,0,1,0) 



3. 



?72 = 



V^^,^ = lim ms 



2P, 



S16 



with 



1 - P'a 



fc'^ = (0, 2VC(i-C), 2C - 1, 



(75) 
(76) 

(77) 
(78) 
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vi = m 



with 



z = 



lim 9^5- 



2P' 



S56 



(79) 



6+6 



(O, v/W, sgn(7}2-77i)(l-27}i)v/C, -2sgn(772~??i)V'7i(l-^i)C) 



5. 



1 

V = - lim V 5y (n2) T, • T^- , 

?2 i, = l 



where 



Pi ■ Pj 



[Pi ■ k)ipj ■ k) ' 

and pi_j is one of pi,p2, 9i, <Z2, fci- are the standard 0] color operators. 



(80) 

(81) 
(82) 



6. 



6 = 6 = A 7}i,2 7^ A f]i^2 7^ 1 A 7)1 ^ 7)2 



In the case of a gluon pair the limit is given by 
V = lim d\s 



(83) 



4 ^ 

^ 2'^*i(|ini)5fei(|2n2)| 



iijkl—l 



Ti}-C^^5,,(Cini, 6^2) 



ij=l 



This approximation is discussed in [Appendix B[ where we also define 5y (fci, ^2). In the case 
of a quark pair we have 



V = Uin 9^5 Tf ^ 1,^(6^1, |2ri2) T, • T, , 

ij = l 



with [57 1 



Ti;j{ki,k2) 



{pi ■ fci) (Pj ■ ^2) + {Pj ■ ki) (pi ■ k2) - {pi ■ Pj) jki ■ k2) 
(fci • k2f [p-, ■ (fci + fcz)] [pj ■ (fci + fcz)] 



(84) 



(85) 



and Pi , in the expressions above is one of pi,p2,qi,q2- 



771 = 6 = A 7)2 7^ A 7/2 7^ 1 A ae = g 



with the collinear parameters specified in Eq. (j76p . whereas p^j- in Eq. 
ki,P2,qi,q2- 



(86) 

is one of pi — 



7)2 = 6 = A 7)1 ^ A 771 7^ 1 A 05 = 5 



Z_r 1 
lim 5H5^^ V -^5,, (77i) T, • Tj 



with the collinear parameters specified in Eq. (j78p . whereas pij- in Eq. 

^^2,^2,91, 92- 



(87) 

is one of pi — 
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m = 572 A ^1 = ^2 = A 7)1,2 7^ A 771,2 1 



S56 



This approximation is discussed in [Appendix C[ and the coUinear parameters are given in 
Eq. dHOl) 

The conditions for the various hmits have been chosen such that only one of the above expressions 
generates a non-vanishing V for a given X. We have also minimized the use of soft limits. For 
example, the double-soft limit contains the double-coUinear and double-soft limit, nevertheless the 
expressions are much lengthier than the double-soft limit of the double-collinear limit, since the 
latter contains no color correlators. This is only a practical choice, since all expressions would give 
the same after taking into account color conservation. 

The double-collinear sector requires the cases 2, 5, 6 and 7 from the triple-collinear sector, as 
well as the following 



1. f/i = A 7)2 = 1 



■xrss s s 
* aia5a2ae 



o pss 2 ^ 

^ lim y\s 

X^O Si5 S26 



(89) 



with the coUinear variables defined in Eq. ((75)) for P^^^^, and in Eq. (1751) for -P^^ae • 



V2 



2 P^^ 
lim fH5^2££ 

X^O S26 



with the coUinear variables defined in Eq. ([75)) . 



?}2 = 1 A 6 = A ?)i 7^ A 7}i 7^ 1 A ag = g 



(90) 



2 

lim d\s^^ 

X^O S26 



4 



11=1 ^1 



m) T, • T 



(91) 



with the coUinear parameters specified in Eq. (I78[) . whereas in Eq. 



is one of pi , P2 " 



Since the V operators are used both for subtraction and integrated subtraction terms, there 
is one more difference of the present approach to the traditional one. There will be transverse 
vectors k± in the integrated subtraction terms. These are usually removed (averaged over) using 
the fact that once in the coUinear limit, one can integrate over the azimuthal angle. Keeping them 
makes our approach less sensitive to simple errors, while not hampering efficiency. 

We do not provide the explicit expressions for all the limits derived according to the rules above. 
On the one hand, it is easy to obtain them, on the other, the formulae are extremely lengthy. We 
believe that it only makes sense to provide a complete, general, working implementation of our 
subtraction scheme. We will return to this in the future. 

There is one more aspect that we can discuss now, namely convergence. Due to the well known 
pointwise nature of the listed limits, when using polarized splitting kernels, the convergence of the 
cross section integrands will be pointwise. We can also assess the rate of convergence. Assume 
that one of the variables x € {t/i, 772, Cii ^2} is rescaled as x ^ k x, k ^ 0, while the others remain 
fixed. If K = implies ^1 = or ^2 or 7)1 = 7)2, but neither 7)1.2 = nor 7)12 = 1, then, ignoring 
logarithmic enhancements. 



(k x) 



(x) 



(92) 
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On the other hand, if k = impUes any of 771 = 0, 7)2 = 0, fji = 1, 7)2 ~ 1, then again up to 
logarithmic enhancements 



(5) 

o 



1 



(k x) « ^ (x) . (93) 



This is the well known inverse square root behavior of collinear limits. The lack of such a behavior 
in the case rji = 7)2 is due to the dependence of the relative angle parameter 773 on the difference 
7)1 ~ fi2, which is quadratic as seen in Eq. (1^ . Due to the iterative derivation of the integrand, 
rescaling several variables leads to a scaling, which can also be obtained iteratively from the above 
formulae. Let us stress, that unless the given limit is a single-soft limit for a final state quark, the 
unsubtracted integrand behaves as 

S^^(k2;) « -E^)(a;) . (94) 

3.2. Leading divergences and leading logs 

We have stressed in the Introduction that one of the main ideas behind our subtraction scheme 
is to avoid any non-trivial analytic integration and perform the entire calculation purely numeri- 
cally. It is, however, advantageous to have at least some analytic formulae to perform tests of the 
implementation. This was already our motivation in deriving the volume of the phase space. It 
turns out that we can also obtain the leading singularity in e directly from the subtraction terms. 
After all, the term corresponds to 771 = 7/2 = = = 0, which means that the reduced 
matrix element is that of the leading order process. The integral in ( is then trivial 

/ dC , =TT. (95) 

Jo VC(W) 

What remains is the two-particle phase space of the leading order cross section. In consequence, 
we obtain 

O (1) , (98) 

O (i) , (99) 

where is the Born cross section for the two channels and can be found in [Appendix D| Later 
on, we will use these formulae to test the normalization and precision of our numerical calculation. 
At this point, we can, however, verify the cancellation of the singularities in the inclusive top 
quark pair production cross section, since the divergences of the other contributions can be found 
in the literature. Indeed, we have 



RR 






ttaa 


RR 






ttaa 


RR 




'^aa- 


>ttqq 


_RR 




qq—>ttq'q' 



2CHC^+4C^)^(^)^,%, + o(^) , (97) 



vv 
"aa^tt 



(4Cl-l-4Ci)l(^)^,%,^-o(l) , (100) 



- + + ■ (101) 

The two color factors in the parentheses in both equations have a different origin. One is given by 
the two-loop virtual corrections and can be read off from the explicit results of 'sgI [s^I in the high 
energy limit (under the assumption that the leading divergence is proportional to the exact Born 
matrix element). Otherwise, it can be obtained from the complete divergence structure presented 
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in [58|. The second color factor is given by the square of the one- loop matrix element, and is easily 
obtained with the help of the I operator from |59| . 

We still need to derive the divergences of the real-virtual corrections. This is much more 
difficult, because the structure of the singular limits of one-loop amplitudes is not the same as 
that of tree-level amplitudes. In other words, it is not enough to just take the term from 
the one- loop amplitude using the I operator and then use the same operator to obtain the 
divergences from the phase space integration. The situation seems to be even more involved, 
because this problem has never been studied for hadronic heavy quark production, where the 
initial states are also partons. Fortunately, the leading singularity is only due to the purely 
massless states and has to factorize, and since it is due to the soft and collinear limit, it is the 
same irrespective of whether we consider initial or final states. In view of these considerations, we 
can use the results obtained for jet cross sections in e^e~ annihilation fl8|. Applying the general 
formulae obtained there, we have 



RV 

-fUto = -2^^HC^ + 8C^)^(g)'af,^,- + o(^) . (103) 
Combining the three contributions, we have indeed 

- O('). (104) 



= 0{-\. (105) 



'qq^tt+X ^ \ 

Proving the cancellation of all divergences is much more difficult, even though we known that 
the cross section is finite by the Kinoshita-Lee-Nauenberg and the factorization theorems (the 
latter to remove initial state collinear divergences). Let us only note, that it is possible to write 
the coefficients of the 1 /e^ singularities through the Born cross sections using convolutions in the 
worst case. Nevertheless, we will refrain from this exercise. 

The knowledge of the leading singularities can also be exploited in another way. Indeed, it 
allows to predict the leading logarithms of (3. It is enough to know the d-dimensional behavior of 
the phase space near threshold, which we have derived in Eq. P7|) as being Z?"-*^"'^. Performing the 
expansion down to the finite part, we obtain 

C-.-.. = i^C^llog^/3 0^.W + C?(log^/3) , (106) 

- ^CHC^ + 4C^)log^/3(^)'4^,- + 0(log^/3) , (107) 

af^^u,, = 0{\oip) , (108) 

= 0(log'/3) • (109) 

Notice that the coefficients are very large. In fact they are about an order of magnitude larger 
than the coefficients of log^ (3 in the total cross section [28] . This will be reflected in the very large 
values of the cross section near threshold. 



4. Implementation 

The implementation of the subtraction scheme described in the previous sections for the par- 
ticular case of top quarks involves a large set of tree level matrix elements, which can, moreover, 
be spin and color correlated. There are many methods to evaluate these in four dimensions. Nev- 
ertheless, we decided to work in conventional dimensional regularization (CDR), thus keeping the 
full (j-dimensional dependence of the matrix elements. A simple way to derive explicit expressions 
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Figure 4: Cut graph classes used for generation of tree level amplitudes squared. The external thick lines represent 
top quarks, the cut thin lines are gluons, ghosts and light quarks, whereas blobs contain trees. Color operators may 
be assigned to each visible line on the left of the cut. Moreover, cut gluon lines may contain spin correlators. 



in such a case is to generate cut graphs for forward scattering amplitudes with up to three loops. 
This procedure based on Cutkosky rules leaves the freedom of the choice of the external states, as 
any configuration can be obtained by crossing. We decided to take the top quarks on the external 
lines of the cut graphs. The tree-level amplitudes may have up to four external gluons, which 
means that there will be a lot of graphs related by symmetry (exchange of the gluon momenta). 
Taking the massless partons as virtual in the cut graphs takes care of the symmetry, in the sense 
that only one configuration is generated and one can symmetrize the complete amplitudes at the 
end. Thus, our expressions are substantially shorter than they would be, if we just squared tree- 
level amplitudes. Our procedure is a simpler version of the old approach of [60], but we would 
prefer the latter if we only had massless partons. The cut graph classes are shown in Fig. |4l We 
insert color operators Tj to all external lines on the left of the cut, with up to four operators in 
total. We take the cut gluon propagators to be in the Feynman gauge and compensate the gauge 
invariance violation by cut ghost lines. We also insert spin correlators according to 



[kl^k\^ + k\^kl^) , (110) 



where the symmetrized version is all that is required in practice, but one may need two different 
transverse vectors as demonstrated in [Appendix A As far as the color correlators are concerned, 



we have exploited color conservation to reduce the number of needed matrix elements. For an 
amplitudes with five partons, we have 

T5IX3) = -(Ti + T2 + T3 + T4)|X3) . (Ill) 

Moreover, since T| = Ca (for all the channels considered) or T| = Cp (in general), we have only 
five correlators to consider 

Ti • T2, Ti • T3, Ti • T4, T2 • T3, T2 • T4 . (112) 

Similarly, in the case of four-parton amplitudes, we have the following two double-correlators 

Ti-T2, Ti-Tg, (113) 

which may be additionally spin correlated. Amplitudes with four partons require, however, also 
quadruple-correlators. At the level of the amplitude 

(Ti • T2)(Ti • T3) 1X2) + (Ti • T3)(Ti • T2) 1X2) , (114) 

nevertheless, since the matrix elements squared are real, we have 

(X2|(Ti •T2)(Ti •T3) 1X2) = (X2|(Ti •T2)(Ti •T3) I7W2)* 

= (X2|(Ti-T3)t(Ti.T2)MA^2) 
= (X2|(Ti •T3)(Ti •T2) IM2) . (115) 
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1,2 




qq — > tt 
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(1,3)(1,3) 
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(1,2)(1,2) 


qq —7- ti 






99 tt 




(1,2)(1,3) 








99 ^ tt 




(1,3)(1,3) 


gq ttq 






99 tt 


1 


(1,2) 


qg ttq 






99 tt 


1 


(1,3) 


qg -J> ttq 







Table 3: The 42 spin and color correlated amplitudes with four, five and six partons, which are needed for the four 
main channels of top quark pair production. 

We thus only evaluate the following correlators 

(Ti • T2)(Ti • T2), (Ti • T2)(Ti • T3), (Ti • T3)(Ti • T3) . (116) 

Notice, that we have not exploited color conservation at the level of unsymmetrized amplitudes. 
This could provide another minor speedup. 

Due to the selector functions and sectors chosen, we are missing phase space. Indeed, following 
Section [231 we have kf > and fcj* > For a given channel specified by flavor assignments to 
the initial and final partons, the additional contributions can be recovered by permuting initial 
and final states independently. Nevertheless, we can use charge conjugation invariance of QCD 
amplitudes together with rotation invariance, which allows to swap pi ^ P2, to reduce the number 
of configurations, which actually need to be evaluated. At the level of six-parton amplitudes, we 
only need the five cases 

99 ttgg, gg ttqq, qq ttgg, qq ttq'q', qq ttq' q' . (117) 

We will not explicitly mention the last amplitude anymore. Interestingly, however, its contribution 
is, for most cases, very close to that with swapped quark and anti-quark. The complete list of 
matrix elements is given in Tab. |3l 

It is interesting to measure the evaluation time of the integrands for a single phase space point. 
The results of such a measurement are shown in Tab.|H The main point of this table is to show that 
the subtraction and integrated subtraction terms are faster in evaluation than the matrix element 
of the double-real radiation process. This is true in all the cases, but the simplest involving six 
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process 


matrix 
element 
[msec] 


using 
Helac 
[msec] 


including 
subtraction 
[msec] 


in quadruple 
precision 
[msec] 


99 tt99 


13 


53 


18 


450 


qq — > ttgg 


0.71 


6.5 


0.81 


27 


99 tM 


0.71 


6.3 


0.97 


32 


qq — >■ ttq'q' 


0.015 


0.52 


0.041 


1.5 



Table 4: Single phase space point timings for the evaluation of the matrix element with or without subtraction. 
The values in the fourth and fifth columns are the worst timings of all the seven sectors. The matrix elements 
of the present implementation have been compiled without optimization, because of the size of the expressions. 
Quadruple precision is obtained using the native implementation of the Intel Fortran compiler. 



quarks. The latter has a very short expression for the amplitude, and it is simply impossible to 
have even simpler subtraction terms. One could, of course, suppose that the relative efficiency 
of the subtraction scheme is simply due to the very inefficient implementation of the six parton 
matrix elements themselves. That this is not the case is shown by comparing our implementation 



with that of Helac [6lll62|. Reversing the argument, one could suppose that Helac is inefficient. 



since our matrix elements have been obtained in the most naive way, and have, moreover, been 
compiled without optimization. This is not true in practice, as Helac uses helicity sampling (or 
rather random polarization vectors), and has not been optimized for spin summed amplitudes, 
which we need here. It would certainly be adva ntag eous to have an implementation of our scheme 
allowing for helicity sampling, as it was done in [63j at NLO. This requires, however, a tremendous 
effort, which we leave for the future. In Tab. SJ we have also quoted timings for computations in 
quadruple precision. In fact, we have used quadruple precision for all the values presented in the 
next section. Due to the cancellations inherent in subtraction schemes, there is always a risk of 
numerical instabilities. In this first study, we have not made any analysis in this direction, and 
decided to avoid the problem altogether using higher precision. This is certainly an issue, which 
requires improvements. We would also like to point out that the implementation of quadruple 
precision of the Intel Fortran compiler, which we used, is rather inefficient and we could have 
gained a speedup factor of at least tree by switching to an external library. 

Let us now discuss our implementation of the phase space integration. The evaluation of the 
cross sections requires seven- and eight-dimensional integrals, corresponding to the two contribu- 
tions from the two-particle phase space Eq. (j46| . The integrands are indeed integrable, but not 
free of singularities. Besides logarithmic singularities due to the e-expansion, there are inverse 
square root singularities as shown in Eq. ((93)) . Our strategy is to use adaptive Monte Carlo inte- 
gration techniques to improve convergence. We do not even bother with remappings for the square 
roots, which could perhaps help, but would be against our approach of avoiding any complications, 
unless, well, unavoidable. Thus, we use Parni [i^] to take care of all singularities. The seven 
sectors are evaluated with the same program in one run, and we even have implemented variance 
optimization d la stratified sampling, but we have not used this feature for the results presented 
here. The main computing time goes into the seven-dimensional integrals, as they involve the 
six-parton amplitudes. The eight-dimensional integrals have smaller contributions and are faster 
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per phase space point. They are thus neghgible as far as resource requirements are concerned. 

There is one more issue connected to numerics that we need to address, and that cannot be 
solved with higher precision. Inherent to subtraction is the fact that the integrands involve a 
cancellation of many digits close to singularities between the matrix elements and their approx- 
imations. If there were no inverse square roots, this might not have been a problem, but high 
precision of the results requires evaluation close to singularities. We thus need to cutoff the phase 
space to avoid instabilities. The form of the cutoff condition can be chosen at will, but since the 
matrix elements grow as l/?7if?2CiC2 7 when the sector variables become small, we have decided to 
require the following 

VI ?72 6 > A . (118) 

The values we chose for A are given in the next section. Here, we evaluate the size the missing 
phase space due to this condition. In fact, one can show that 

<5<i>„(A) = f[ da, (j[ a, < A j = A ^ log' (^1) . (119) 

In our case n — A. The size of the missing phase space cannot be used to estimate the implied 
error on the total cross section, due to the presence of square roots and logarithms. In practice, 
one evaluates the cross section at two values of the cutoff, say A and A', with A' < A. It is 
expected that the improvement in the error will be of the order of A/A'. 

The final issue we would like to discuss is that of testing. In a project this size, it is important 
to check correctness at all stages. We have performed the following checks 



1. 



Phase space volume We have evaluated the volume analytically in Section [2?T1 By setting 
the matrix elements to unity, we can obtain a numeric result with our implementation. For 
example, with 

s = 2, m=l, ^1 = 1, (120) 



the exact value reads 



3e 



(^5.9719120 + 87.151375 6 + 630.84755 
+3019.2212 + 10746.200 e"*) x 10"^^ , (121) 
whereas 10 000 000 generated points suffice to obtain 

J d<i>4 w ((5.9697 ± 0.0028) + (87.131 + 0.033) e+ (630.77+ 0.18) 



47r 



^(3019.2 + 0.64) 6^ + (10747+ 2.1) e-^) x 10"^^ . (122) 



e-contribution to the two-particle phase space Section [231 contains a result for the integral 

2 



of (fci • (7i)^. This integral is the simplest object sensitive to ^$2'''. Unfortunately, the 



contribution is tiny at best. The analytic result for s ~ 1,to^ = /.t^ = 2/9 is 

-J^ r / d$4 1 ^— ^ I = f 1.4553533 + 10.106976 e + 34.956382 

-P4(s, e) J \ s ) V 

+80.126733 e- + 136.49975 A. x 10"^ , (123) 
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whereas with very high statistics, we have 



P4(s,e) 



^^{d\e) ^ j ^ (^(1.45536 ± 0.000014) + (10.1056 ± 0.000085) e 

+ (34.9482 ± 0.00035) + (80.1058 ± 0.0017) 
+ (136.473 ± 0.0084) e"*) x 10^^ , (124) 



/ f j = ((0.001413 ± 0.000011) e + (0.008145 ± 0.000056) 



P4{s,e) 



+ (0.02064 ± 0.00013) + (0.02469 ± 0.00020) j x 10"^ . 

(125) 

Clearly, there is disagreement between the ^$4'^''^'' values and Eq. (jl23l) starting from order 
e. Besides the highest order, e^, the difference is at least at the level of lOcr. Together with 
the e-contribution we obtain, however, 



Pi{s,e) 



J d^4 (^'^^^-^^ = (^(1.45536 ±0.000014) + (10.1070+ 0.000086) e 

+ (34.9564 ± 0.00035) + (80.1265 ± 0.0017) 
+ (136.498 ± 0.0084) e"*) x 10"^ . (126) 



3. Leading order cross sections With a jet function that forces all massless partons to be sep- 
arated, we obtain the leading order cross section for top quark pair production in association 
with two jets. We use the following simple setup 





.2 ■ Pi, 2 


> 10""* s , fci • fca > 10~^ s , 


(127) 


together with 








Ecu = 400 GeV 




172.6 GeV , a,(mt) = 0.107639510785815 . 


(128) 


We obtain 










^ttgg 


= (7.75 + 0.018) X 10^2 j^b , 


(129) 


^gg- 


-^ttqq 


= (4.81 + 0.019) X 10""* nb , 


(130) 


CTqq- 


-^tigg 


= (2.86 + 0.0065) X 10^2 j^b , 


(131) 


CTqq^ 


tiq' q' 


= (3.55 + 0.010) X 10""* nb , 


(132) 


whereas Helac gives 








'^gg- 


^tigg 


= (7.76 + 0.039) X IQ-^ nb , 


(133) 


'^gg- 


-^tiqq 


= (4.77 + 0.025) X lO^Sib , 


(134) 


<Jqq- 


-^tigg 


= (2.88 + 0.0078) X 10^2 nb , 


(135) 


^qq-> 


tiq' q' 


= (3.54 + 0.019) X 10""* nb . 


(136) 



We point out that the results have been obtained with 1 000 000 Monte Carlo events, and 
our phase space parameterization seems to be more efficient than that of Helac. 
We have also checked our matrix elements for the six-parton processes at chosen phase space 
points against those obtained with Helac. We have reached agreement within expected 
numerical accuracy. 
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Pointwise convergence in d-dimensions We can test the scaling given in Section 13.11 in all 

the limits, while expanding the results up to e*^. This has the advantage of showing that 
the subtraction terms are correct to all orders in e, since the highest order of expansion of 

the matrix elements is e**. The further terms of the expansion of follow from the 

e-dependence of the measure. We have evaluated all the fifteen limits of the seven sectors 
of the four main processes, which amounts to 420 cases. Needless to say, we obtained the 
correct behavior. To illustrate the tests, we only give one case, which sector of the process 
gg — >■ ttgg for the following configuration 



1 , f3 



C 



1 

3 ' 



m 



1 



1 



cos = g > 



X 10 



1 



-40 



6 



1 



X 10 



-40 



- , cos PQ = 1 



(137) 



The results quoted correspond to = 1, /i = m, and we did not include the factors for spin 
and color averages, identical final states and fiux. Without subtraction the integrand is 



2.92 X 10^^ + 1.14003 x 10"" e 



+2.89548 X 10*'^ 



,82 

2.82539 X 10 



2.22528 X 10^'' 



88 g4 



(138) 



The size of the coefficients is due to the expected growth of the integrand proportional to 
the inverse of the product of the four sector variables, which in this case is of the order lO®'^. 
Once subtraction and integrated subtraction terms are included, the integrand becomes 



9.98769 X 10" 



1 



0.0000869195 



1 



1 



1 



-40.9821 + 1961.23 e + 75083.4 + 2.39952 x 10^ e 



0.00676374 — -f 0.641465 

3 ' 6.61228 X 10^ e"* 



(139) 



We observe no growth of the coefficients, since the rescaling of 772 and ^2 corresponds to 



, there is one that cancels the 



the limit 771 ~ 7)2 and ^^2 = 0. Amongst all terms in 

divergence of E'*^!^ alone, and contains the limit of Wig at 771 = 7)2 and ^2 = 0. We can 
compare it order by order in e to S^'^*). The result is 



.{Si) 



1 - 



-'approx I 



= -2.68 X 10" 



2.68 X lO"''^ e 



2.68 X lO^'^^ 



-2.67 X 10"*^ - 2.67 x lO"**^ 



(140) 



The very high numerical precision of these tests is important, because minor mistakes can 
sometimes show, for example, at the eight digit and would be impossible to find with a 
double precision implementation in Fortran. The above numbers have been obtained with 
an arbitrary precision implementation in Mathematica. We have also compared quadruple 
precision results from Fortran with those of Mathematica at ordinary points making sure 
that both implementations agree. 

Further tests, most notably the agreement of the numerical values for the leading divergence with 
the analytic formulae from Section 13.21 as well as cutoff independence are described in the next 
section. 

Finally, let us point out that this project has only been possible thanks to the use of complicated 
external software systems. These are listed in [Appendix E[ 
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Table 5: Coefficients of the Laurent expansion of the fqq^ttgg function. 

5. Results 

We are now ready to present our numerical results for the Laurent expansions of the cross 
sections. We turn to dimensionless functions of the velocity /3 with the help of the fohowing 
definition 

4 

'^ab'lticdis^ = m^.as, e) = e) , (141) 

where a, b are initial, and c, d final state massless partons. As shown on the left hand side above, 
we only give values for the case fi — m, since the dependence on the scale can be obtained from 
renormalization group equations. The / functions admit the following expansions 



fat^tUl3, e)^Yl ''faLuM + 0{e) ■ (142) 
We sample the functions at twenty equidistant points within the /3 variation range [0, 1] 

/3, = ^^, *-l,...,20. (143) 

This is most probably sufficient to obtain a decent fitting function as has been done in the classic 
paper [g^. To this, we add one point very close to threshold, /3 = 0.001, and one very close to the 
infinite energy limit /3 = 0.999. As we will demonstrate, the study of the latter limit will require a 
denser sampling in the case of gluonic initial states, even though this is not immediately relevant 
to phenomenology due to the limited /3 range of current colliders producing top quarks. 
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Table 6: CoefScicnts of the Laurent expansion of the fqq^ttgg function. 

As we are interested in total cross sections at first, we use a trivial jet function equal to one. 
For our final results, we use a cutoff of A = 10^^. This implies a missing volume of 

J$4(10"^) = 8.4 X 10^^ . (144) 

In order to test the independence of the results from the cutoff, we will also use a higher value of 
A = 10~^, which amounts to 

J$4(10"^) = 5.5 X 10^4 _ (-145) 

These numbers can only be used as estimates of the relative integration errors implied by the 
cutoff in the non-singular case as discussed in Section |4l Nevertheless, the improvement obtained 
by lowering the cutoff by an order of magnitude should amount to a factor of three, which will 
allow us to provide realistic estimates of the quality of the final results. In fact, we will aim at a 
situation, in which we will be restricted by the integration error rather than by the cutoff. 

For each cross section and each value of /3, our simulations are performed with a total sample 
of 10 000 000 generated Monte Carlo events. This number is very close to the number of accepted 
events, since the only restriction on the phase space is given by the small cutoffs. The quality 
of the obtained results with this sample is discussed in the following. Nevertheless, we justify 
its fixed size by the purpose of this publication, which is to prove the usefulness of our approach 
and provide relevant numbers, while not necessarily giving the highest quality estimates. The 
latter exercise is left for a future publication containing complete cross sections, and not only the 
double-real contribution. We stress already here, that there is always a risk of underestimated 
errors with complicated integrands and low relative errors. Therefore, our results can probably 
only be trusted up to an additional factor multiplying the quoted errors. A safe bet with no special 
justification would be a factor of two. This is unavoidable, as practice shows, and the only way 
not to have to worry about the errors is to substantially increase the statistics to the point, where 
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the estimated precision will vastly surpass the practical requirements. This will of course be done 
in the mentioned future publications. 

Finally, we re-stress that all computations have been performed in quadruple precision, in order 
to remove at least one source of concern, namely numerical instabilities. Using higher precision is 
by itself not yet reassuring enough, but two calculations with two different cutoffs and agreement 
within expectations should be. 

We start our presentation with the qq — ^ ttgg channel, as it contains all the complications as 
far as the singularity structure is concerned, but can be evaluated in a shorter time in comparison 
with the cutting edge gg — ?> tigg process, due to the faster per phase space point computation. 
The results for the expansion coefficients can be found in Tabs. [5] and [6l Notice that we have 
not given the values for the leading singularity as it is known analytically from Section 13.21 
Striking are, of course, the very large values close to threshold, which were, however, expected 
from the leading logarithmic behavior determined in the very same Section 13.21 At this point, 
we can comment on the integration errors. It is interesting that, while the relative error varies 
substantially, the absolute error stays more or less the same up to some factor. The reason for this 
behavior that we will see in all subsequent results is that the dominant, logarithmic contribution 
in /? has much higher precision, due to simpler functional dependence (neither logarithms nor 
inverse square roots in the integration parameters). The remaining functional dependence is due 
to cancellations between contributions of the different sectors, which in turn have all more or less 
a similar absolute error. If we now assume that the error is a constant 2 x 10~^ (not exactly 
an upper bound, but rather a realistic estimator for integration with the partonic flux), we can 
obtain the implied relative error on the top quark pair production cross section at the TeVatron, 
where this channel dominates. It turns out, that the error would be 0.2%, which is more than 
acceptable. 

Let us now study the cutoff dependence of the results. The latter is illustrated in Tab. [TJ 
where we give the difference between the finite parts of the cross section evaluated at both cutoffs, 
10"'' and 10~^, compared to the sum of the integration errors. For many values of /3, we notice 
that the difference is larger than the errors, albeit not by a large factor. Taking into account that 
we can consider the difference to be close to the actual variation with the cutoff, when changing 
from 10"^ to 10"^, and that we expect a variation smaller by a factor of three, when stepping 
to 10~®, we expect that the actual error due to finite cutoff is lower than the integration error, 
as certified by dividing all numbers by three and comparing to Tab. |6l The table also contains 
the e-contribution to the two-particle phase space. We note, that it is at the level of the current 
integration error for most points. 
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Figure 5: Difference between the value of the coefficient of the leading divergence of fqq^ttgg obtained by numerical 
integration and the exact expression Eq. I I97II normalized to the value of the latter. On the left panel, the lower 
cutoff A = fO~^ has been used, whereas on the right, A = 10~®. Notice the shifted scale on the vertical axis 
between the left and right panels. 
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Table 7: A comparison between the finite parts of the fqq^ttgg function evaluated at two different values of the 
cutoff A defined in Eq. I|118|l . A = 10~^ and A = 10~®. The second column represents the difference between the 
values, whereas the third the sum of the integration errors. The fourth column contains the e-contribution of the 
two-particle phase space. 



In order to check the normahzation of our results, we can also compare the numerical estimate 
of the leading singularity, with the prediction from Section [X^ This is done in Fig. [S]for 

both values of the cutoff, where we plot 

I (47r)2/<rtl- -2CF(CA+4Ci.)/*°' - 

where fj^^l^^i is defined by the Born cross section in [Appendix D[ We notice that in the case of 
the higher cutoff, A = 10~^, the cutoff dependence is noticeable due to the tiny numerical errors. 
What is slightly more worrisome is that a few errors are indeed underestimated. In this case, this 
is due to the fact that the integration errors are very small, below permille, but the function has 
integrable singularities in ( at the integration boundaries, and is thus not entirely well behaved. 
Based on this, we can only expect worse from the much more singular finite parts of the cross 
section. 

Finally, in Fig. [SI we show the cross section after subtracting the leading logarithm in /3. The 
latter would make the plot span an order of magnitude more. Notice that on the scale of this plot, 
the integration errors would not be noticeable. 

We can, in principle, repeat the same discussion for the most complicated and computationally 
intensive channel, gg — >■ ttgg. The numerical values are given in Tabs. [S] and [31 where we have 
again omitted the leading singularity. The essential difference to the quark annihilation channel is 
in the about twice larger absolute errors of the finite part for most of the /3 range. One can again 
estimate that if the error is consider constant and equal to 4 x 10^2^ then the implied uncertainty 
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Figure 6: Finite part of fqq^ttgg after removing the dominant logarithmic term Eq. I I107I I. 

at the LHC would be 0.4%, which is also more than acceptable. As before, we can demonstrate 
that the cutoff dependence is lower than the numerical integration error. 
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Figure 7: Finite part of fgg^ttgg after removing the dominant logarithmic term Eq. 11061 1. 

In Fig. [71 we show the finite part of the cross section after removing the dominant logarithm 
in /3. We notice the very steep rise at the end of the range. This phenomenon is well known from 
the next-to-leading order cross section in the gluon fusion channel. For our numerics it has the 
unpleasant feature of making the calculation slightly unstable due to the extreme sensitivity to 
the value of /3. This also points to an underestimated error a.t /S — 0.999, anyway quoted to be 
rather large in Tab. [SI None of this is relevant to phenomenology at present, but we will try to 
get a better handle of this problem in the future. 

The remaining two cross sections are much less interesting. They are less singular both in e 
and log /3, and are moreover much smaller even after multiplication by the number of massless 
quark species. We give the numbers in Tabs. [TOl [HI [El and [131 and show the respective finite 
parts in Fig. \8\ 
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Table 8: Coefficients of the Laurent expansion of the fgg^ttgg function. 

6. Conclusions 

The purpose of this work was to prove the usefulness of the Stripper approach to the problem 
of double-real radiation. We have considered the phenomenologically relevant case of top quark 
pair production, and evaluated the cross sections for the dominant channels. We have given most 
of the formulae needed for the implementation, and demonstrated pointwise convergence and 
efficiency. The immediate consequence is the possibility to evaluate the complete cross sections 
for top quark pair production after inclusion of the double- virtual, and real- virtual contributions. 
Although this requires quite some effort, we do not see any conceptual problems, unlike in the 
present case of double-real radiation. 

We would like to point out that there are three directions of further development. First, there 
are many technical improvements of our implementation that can be studied. The most important 
are the analysis of numerical instabilities and implementation of a more efficient multi-precision 
library, the latter being almost trivial. Although not absolutely necessary, this work is always 
part of software maturation in the case of higher order calculations. The second direction involves 
applications to similar process, in particular removing some of the subtraction terms is sufficient 
to treat e~^e~ — )■ tt + X, and pp{p) — W~^W~ +X (and other gauge boson final states) at NNLO. 
The last, and probably the most interesting, direction is application of Stripper to final states 
with massless particles, such as dijet production. This requires the specification of the phase space 
in the case of initial and final state singularities, but as noticed in [l| involves the same treatment 
of the unresolved partons. 
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Table 9: Coefficients of the Laurent expansion of the fgg^ttgg function. 
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Appendix A. Collinear limits and splitting functions 

In this appendix, we will reproduce the sphtting functions that have been used in the derivation 
of the subtraction terms. The formulae are taken literally from [s^l (see also 6^ 67[). We start 
by defining the notation for the matrix elements 



■' 'ai,a2,... 



(A.1) 



where the Si indices stand for spin, the Ci for color, and the for parton flavor. With this object, 
we define spin correlated amplitudes squared 



spins T^si ,s'^ colors 



C2,...;Si,S2, 
,02,... 



■{pi,P2,...) M'a\'^a2.:f''''^''''{Pl,P2, 



t 



(A.2) 

Having defined the matrix elements, we now turn to next-to-leading order collinear limits of 
amplitudes. We first define the limits through auxiliary vectors 



Pi = zp^ + /Cj 



Sl2 =2pi-p2 = - 



z 2p ■ n 

z{l-z) 



, p^ = (l-zK-fc'[ 



fc_L ^ , 



1 ~ z2p ■ n 



(A3) 
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± 


8.6 X 


10 


6 


0.675 


-1.208 X 


10 


3 


± 


1.2 X 


10 


6 


-6.310 X 


10 


3 


± 


8.3 X 


10 


6 


0.725 


-1.256 X 


10 


3 


± 


1.3 X 


10 


6 


-4.836 X 


10 


3 


± 


8.0 X 


10 


6 


0.775 


-1.264 X 


10 


3 


± 


1.3 X 


10 


6 


-2.797 X 


10 


3 


± 


7.6 X 


10 


6 


0.825 


-1.222 X 


10 


3 


± 


1.3 X 


10 


6 


-8.582 X 


10 


5 


± 


7.4 X 


10 


6 


0.875 


-1.107 X 


10 


3 


± 


1.2 X 


10 


6 


+3.444 X 


10 


3 


± 


7.9 X 


10 


6 


0.925 


-8.799 X 


10 


4 


± 


9.5 X 


10 


7 


+8.174 X 


10 


3 


± 


9.7 X 


10 


6 


0.975 


-4.551 X 


10 


4 


± 


6.3 X 


10 


7 


+1.546 X 


10 


2 


± 


1.3 X 


10 


5 


0.999 


-3.921 X 


10 


5 


± 


5.9 X 


10 


7 


+2.366 X 


10 


2 


± 


1.9 X 


10 


5 



Table 10: Coefficients of the Laurent expansion of the fgg-ntqq function. 

where p'^ = n'^ = p ■ k± — n ■ k± = 0. Notice that all vectors here and below are outgoing. The 
case we are interested in, namely some of the vectors being in-going, is recovered by crossing. In 
the above collinear limit, the matrix element factorizes as follows 



Sl2 

The splitting functions Pafa, depend on the parton flavors. For the general case 



they read 



a{p) ai{zp + k± + 0{kl)) + a2((l - z)p -k± + C>(fci)) , 



{z, kr,e) = P?^ {z, kr,e) = Sss' Cf 



1 + z^ 
1-z 



-e{l-z) 



(A.4) 



(A.5) 



(A.6) 



P^^ {z, kr,e) = P-^l (z, kr.e)= 5ss' Cf 
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1 + (1 - zf 



(A.7) 



P^^{z, k±;e)= P^^{z, k±;e)= Tf 
Pi;^^{z,kr,e) = 2CA 



-g^^+Az{l-zf-^ 



-r ( + V) - 2(1 - e).(l - z) 



\-z 



k^, 



(A.8) 
(A.9) 
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/3 






0.001 


-3.622 


X 


10-^ 


± 


6.2 


X 


10-" 


-8.500 


X 


10 


2 


± 


1.8 


X 


10-* 


0.025 


-2.672 


X 


10-2 


± 


1.9 


X 


10-5 


-3.366 


X 


10- 


1 


± 


2.7 


X 


10-4 


0.075 


-4.108 


X 


10-2 


± 


3.2 


X 


10-5 


-3.636 


X 


10- 


1 


± 


3.4 


X 


10-^ 


0.125 


-4.604 


X 


10-2 


± 


4.1 


X 


10-5 


-3.275 


X 


10- 


1 


± 


4.4 


X 


10-* 


0.175 


-4.781 


X 


10-2 


± 


4.3 


X 


10-5 


-2.869 


X 


10- 


1 


± 


3.7 


X 


10-* 


0.225 


-4.806 


X 


10-2 


± 


4.6 


X 


10-5 


-2.495 


X 


10- 


1 


± 


3.7 


X 


10-^ 


0.275 


-4.742 


X 


10-2 


± 


4.8 


X 


10-5 


-2.173 


X 


10- 


1 


± 


3.5 


X 


10-* 


0.325 


-4.607 


X 


10-2 


± 


5.0 


X 


10-5 


-1.881 


X 


10- 


1 


± 


3.6 


X 


10-* 


0.375 


-4.404 


X 


10-2 


± 


5.2 


X 


10-5 


-1.612 


X 


10- 


1 


± 


3.5 


X 


10-^ 


0.425 


-4.137 


X 


10-2 


± 


5.5 


X 


10-5 


-1.364 


X 


10- 


1 


± 


3.6 


X 


10-4 


0.475 


-3.770 


X 


10-2 


± 


5.3 


X 


10-5 


-1.111 


X 


10- 


1 


± 


3.5 


X 


10-4 


0.525 


-3.296 


X 


10-2 


± 


5.7 


X 


10-5 


-8.534 


X 


10- 


2 


± 


3.8 


X 


10-4 


0.575 


-2.716 


X 


10-2 


± 


5.9 


X 


10-5 


-5.918 


X 


10- 


2 


± 


3.7 


X 


10-4 


0.625 


-2.032 


X 


10-2 


± 


5.8 


X 


10-5 


-3.295 


X 


10- 


2 


± 


4.0 


X 


10-4 


0.675 


-1.234 


X 


10-2 


± 


5.7 


X 


10-5 


-7.767 


X 


10- 


3 


± 


3.7 


X 


10-4 


0.725 


-3.706 


X 


10-3 


± 


5.2 


X 


10-5 


+ 1.446 


X 


10- 


2 


± 


3.3 


X 


10-4 


0.775 


+5.117 


X 


10-3 


± 


5.0 


X 


10-5 


+3.236 


X 


10- 


2 


± 


3.1 


X 


10-4 


0.825 


+1.331 


X 


10-2 


± 


5.6 


X 


10-5 


+4.375 


X 


10- 


2 


± 


2.9 


X 


10-4 


0.875 


+1.891 


X 


10-2 


± 


7.8 


X 


10-5 


+4.700 


X 


10- 


2 


± 


2.9 


X 


10-4 


0.925 


+1.833 


X 


10-2 


± 


4.4 


X 


10-5 


+4.423 


X 


10- 


2 


± 


2.4 


X 


10-4 


0.975 


-3.231 


X 


10-3 


± 


6.0 


X 


10-5 


+6.853 


X 


10- 


2 


± 


2.6 


X 


10-4 


0.999 


-6.995 


X 


10-2 


± 


1.9 


X 


10-^ 


+3.312 


X 


10- 


1 


± 


1.2 


X 


10-3 



Table 11: Coefficients of the Laurent expansion of the fgg-ntqq function. 



Let us now turn to the more complicated case of triple-collinear limits. Consider the set of 
three vectors 

1,2 „^ 

P^ = XiP>' + e^i-^^ , i = 1,2,3, (A.IO) 

where as before = n? = p - k±i = n ■ k±i = 0. This configuration fulfills no other constraints, 
but rather the limits are expressed through derived variables 

Zi = ^,3^ , (A.11) 
Z^j=i Xj 

= k^i--^j2K, . (A.12) 



We also define 



J. _ o ZjSjk ZjSjk Zi Zj ( \ TQ\ 

Zi ~r Zj z% ~r Zj 



with Sij = {pi +}3j)2. 

The factorization formula is now 

|Mai.a.,a3,...(Pl,P2,P3,...)P^ f^^V-^) V^ixP, ■ ■ ■) P^a.a, , (A.14) 

with si23 = {pi +P2+ Ps)^ and x = xi + X2 + xs- 
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/3 


e-3 


e-2 


0.001 


-1 


962 


X 


lo-*^ 


± 


8.8 


X 


10-4U 


-1.418 


X 


10-4 


± 


6.5 


X 


io-« 


0.025 


-4 


907 


X 


10-5 


± 


2.4 


X 


10-8 


-1.967 


X 


10-3 


± 


9.6 


X 


10-7 


0.075 


-1 


464 


X 


10-4 


± 


1.4 


X 


10-7 


-4.255 


X 


10-3 


± 


2.5 


X 


10-^ 


0.125 


-2 


406 


X 


10-4 


± 


1.4 


X 


10-7 


-5.749 


X 


10-3 


± 


3.1 


X 


10-6 


0.175 


-3 


297 


X 


10-4 


± 


1.9 


X 


10-7 


-6.740 


X 


10-3 


± 


3.7 


X 


10-6 


0.225 


-4 


129 


X 


10-4 


± 


2.5 


X 


10-7 


-7.342 


X 


10-3 


± 


4.2 


X 


10-6 


0.275 


-4 


866 


X 


10-4 


± 


3.0 


X 


10-7 


-7.593 


X 


10-3 


± 


4.5 


X 


10-6 


0.325 


-5 


508 


X 


10-4 


± 


3.4 


X 


10-7 


-7.558 


X 


10-3 


± 


4.7 


X 


10-6 


0.375 


-6 


029 


X 


10-4 


± 


3.9 


X 


10-7 


-7.250 


X 


10-3 


± 


4.7 


X 


10-6 


0.425 


-6 


429 


X 


10-4 


± 


4.2 


X 


10-7 


-6.723 


X 


10-3 


± 


4.7 


X 


10-6 


0.475 


-6 


682 


X 


10-4 


± 


4.5 


X 


10-7 


-5.998 


X 


10-3 


± 


4.3 


X 


10-6 


0.525 


-6 


783 


X 


10-4 


± 


4.6 


X 


10-7 


-5.117 


X 


10-3 


± 


4.0 


X 


10-6 


0.575 


-6 


724 


X 


10-4 


± 


4.5 


X 


10-7 


-4.113 


X 


10-3 


± 


3.5 


X 


10-6 


0.625 


-6 


503 


X 


10-4 


± 


4.4 


X 


10-7 


-3.037 


X 


10-3 


± 


3.1 


X 


10-6 


0.675 


-6 


124 


X 


10-4 


± 


4.2 


X 


10-7 


-1.945 


X 


10-3 


± 


2.6 


X 


10-6 


0.725 


-5 


569 


X 


10-4 


± 


3.9 


X 


10-7 


-8.748 


X 


10-4 


± 


2.2 


X 


10-6 


0.775 


-4 


869 


X 


10-4 


± 


3.6 


X 


10-7 


+9.961 


X 


10-5 


± 


2.0 


X 


10-6 


0.825 


-4 


001 


X 


10-4 


± 


3.1 


X 


10-7 


+9.114 


X 


10-4 


± 


1.7 


X 


10-6 


0.875 


-2 


998 


X 


10-4 


± 


2.4 


X 


10-7 


+ 1.461 


X 


10-3 


± 


1.5 


X 


10-6 


0.925 


-1 


874 


X 


10-4 


± 


1.5 


X 


10-7 


+ 1.609 


X 


10-3 


± 


1.3 


X 


10-6 


0.975 


-6 


459 


X 


10-5 


± 


4.3 


X 


io-« 


+ 1.052 


X 


10-3 


± 


6.0 


X 


10-7 


0.999 


-2 


617 


X 


10-6 


± 


2.0 


X 


10-9 


+ 1.013 


X 


10-4 


± 


1.1 


X 


10-7 



Table 12: Coefficients of tfie Laurent expansion of the fqq^ttq'q' function. 



The complete set of splitting functions is (in the case of spin conservation, we give only the 
spin averaged splitting functions (^^10203)) 



I1 92 93 / 



1 ^ „ S123 

2 S12 



^12,3 , 4Z3 + (Zi - Z2f 



S12S123 



Zl + Z2 



(1 -2e)\z, + Z2- — 
S123 



(A.15) 



Notice that we have omitted the case of identical quarks, which is not needed in the present paper. 
The splitting functions for this case can be found in |57| . The remaining functions are 



919293' ' 



(A.16) 



with 



/ p(ab) \ ^ 
\ .91.9293/ 



'123 



I 2si3S23 



23 



1 + Z 



2 

3 _ g^l 



2? + zl 



Z1Z2 



Z\Z2 



-e(l + e) 



5123 
Sl3 



Z3(l - Zl) + (1 - Z2)3 , 2n , ^ ^ 2 , ,2^1-^2 
h e (1 + Z3) - t{z-^ + Z1Z2 + Zj)' 



2122 



X S23 
■513. 



Zl22 



+ (1^2) , 



(A.17) 



36 



/3 




eO 


0.001 


-5.122 


X 


10-3 


± 


2.4 


X 


10-" 


-1.232 


X 


10-^ 


± 


6.2 


X 


10-5 


0.025 


-3.933 


X 


10-2 


± 


2.0 


X 


10-5 


-5.229 


X 


10-1 


± 


2.9 


X 


10-1 


0.075 


-6.148 


X 


10-2 


± 


3.6 


X 


10-5 


-5.894 


X 


10-1 


± 


6.4 


X 


10-4 


0.125 


-6.810 


X 


10-2 


± 


3.9 


X 


10-5 


-5.340 


X 


10-1 


± 


4.0 


X 


10-4 


0.175 


-6.802 


X 


10-2 


± 


4.1 


X 


10-5 


-4.529 


X 


10-1 


± 


3.6 


X 


10-4 


0.225 


-6.408 


X 


10-2 


± 


4.1 


X 


10-5 


-3.670 


X 


10-1 


± 


3.5 


X 


10-4 


0.275 


-5.768 


X 


10-2 


± 


4.3 


X 


10-5 


-2.848 


X 


10-1 


± 


3.8 


X 


10-4 


0.325 


-4.992 


X 


10-2 


± 


3.7 


X 


10-5 


-2.112 


X 


10-1 


± 


2.6 


X 


10-4 


0.375 


-4.117 


X 


10-2 


± 


3.7 


X 


10-5 


-1.453 


X 


10-1 


± 


2.4 


X 


10-4 


0.425 


-3.225 


X 


10-2 


± 


3.8 


X 


10-5 


-9.050 


X 


10-2 


± 


3.4 


X 


10-4 


0.475 


-2.344 


X 


10-2 


± 


2.7 


X 


10-5 


-4.580 


X 


10-2 


± 


1.7 


X 


10-4 


0.525 


-1.523 


X 


10-2 


± 


2.9 


X 


10-5 


-1.203 


X 


10-2 


± 


2.9 


X 


10-4 


0.575 


-7.844 


X 


10-3 


± 


2.6 


X 


10-5 


+ 1.230 


X 


10-2 


± 


1.7 


X 


10-4 


0.625 


-1.715 


X 


10-3 


± 


2.1 


X 


10-5 


+2.634 


X 


10-2 


± 


1.6 


X 


10-4 


0.675 


+2.967 


X 


10-3 


± 


2.0 


X 


10-5 


+3.260 


X 


10-2 


± 


1.1 


X 


10-4 


0.725 


+5.959 


X 


10-3 


± 


1.4 


X 


10-5 


+3.152 


X 


10-2 


± 


8.3 


X 


10-5 


0.775 


+7.193 


X 


10-3 


± 


1.3 


X 


10-5 


+2.558 


X 


10-2 


± 


9.4 


X 


10-5 


0.825 


+6.576 


X 


10-3 


± 


1.2 


X 


10-5 


+1.680 


X 


10-2 


± 


6.1 


X 


10-5 


0.875 


+4.164 


X 


10-3 


± 


7.8 


X 


10-*^ 


+8.379 


X 


10-3 


± 


4.3 


X 


10-5 


0.925 


+3.782 


X 


10-4 


± 


8.4 


X 


10-'^ 


+3.761 


X 


10-3 


± 


4.1 


X 


10-4 


0.975 


-3.322 


X 


10-3 


± 


2.9 


X 


10-^ 


+9.508 


X 


10-3 


± 


1.3 


X 


10-5 


0.999 


-1.074 


X 


10-3 


± 


8.8 


X 


10-7 


+7.280 


X 


10-3 


± 


7.3 


X 


10-^ 



Table 13: Coefficients of the Laurent expansion of the fqq^ttq'q' function. 



/p(nab) \ 

\ (7K72 03/ 



Similarly 
with 



+ 



(1-e) 



■-12,3 



+ 



4s?2 ' 4 2 



+ 



'123 



2S12S13 



(l-Z3)'(l-e) + 2z3 

Z2 



z|(l-e) + 2(l-Z2) 



1-^3 



'123 



4S13S23 



Z3 



(l-Z3)'(l-e) + 223 



2:1-22 



+ e(l-e) 



5123 
2S12 

S123 

2^ 



(1-e) 



23(-I 



_ (1 - Z2f + Zl-Z2 _ 
^ ' 22(1-23) 



22(1 - 23) Z2{1 - 23) 

2(1 - 22)(22 - 2:3) 



22(1 - 23) 



23(1 - 21) + (1 - Z2Y 
21 22 



+ e(l-22) 



2|_+2| 
21 22 



■21+22 



(lo2) . 



5192^3 ^^-^-^-^519293 ^^ -^ 5192^3 ' 



^fiv (nab) 



(A.I8) 



(A.19) 



P. 



t^iy (ab) 



+ 



45123 
S12S13 



^ ^ 2si23S23 + (1 - e)(si23 - S23)' 



S12S13 



^^3^2 "I" ^2^3 ~ 0- ~ ^)^l'^l^ 



(A.20) 
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Figure 8: Finite parts of fgg^ttqq (left) and fqq^ttq'q' (right). 



A/^z^ (nab) 
3192 93 



Finally 



pfiu 

913293 



1 I Sl23 



'23 



Si 23 



- 16- 



'2^3 







( k2 




\Z2 




\Z2 " 





Si 23 
S12S13 



251233^'' - 4(fc2^fc3^ + kt;k!^ - (1 - e)A:f fcr) 



Sl23 
S12S23 



(1- 2e) + 2 
2S123.9'' 



Sl23 1 - Z3 

si2 2;i(l - Zi[ 
.^2(1 -2zi) 



S123 1 - zi + 2z1 



S23 



^1(1 -zi) 

.v2 



Zl(l - Zl) 



Zl(l - Zl) 



+ 8{1 -e)k^k^ 



mk^ + e,k^,) (^ll^^T^ + (1 - ^)) 



+ (2 o 3) 



(A.21) 



4s?2 



-5'''i?2,3 + 16s 



2^2 



123- 



1^2 



/fc2 


4)' 


^ fc2 




U2 




1^2 


Zl 1 



4 S12 Z3 



2(1 - 



4z| 



1 



Sl23(l - e) 

S12S13 

S123 
2(1 -e)-' 



I /bo k' 



r. 1-2^3 



2:3 



1 - 2z3(l - Z3) 

Zi(l - Zl) 
u 1-2Z2 



^2 Z3(l-Z3)^'^=^'^^ Z2(l-Z2; 

/^4z2Z3 + 2zi(l-zi)- 1 l-2zi(l-zi) 



(1-Z2)(1-Z3) 
/2Z2(1-Z2) 



V 2 3 3 2y ^ Z3(l-Z3) 



Z2Z3 

(5 permutations) 



(A.22) 



Appendix B. Soft limits in the presence of massive partons 



While Ref. 57| contains a summary of the behavior of QCD matrix elements in singular limits at 
next-to-next-to-leading order, the authors restricted themselves to the case of massless partons. 
Since massive partons do not induce coUinear singularities, we need only consider the soft limit. 
It is well known that the eikonal current has the same form in both massless and massive cases. 
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This implies that as long as we describe strongly ordered limits, no modification of the expressions 
is needed. Surprisingly, one observes a difference in the double-soft limit, in which the energies of 
both gluons (there is nothing special in the case of a soft quark pair) vanish at the same rate. To 
be more specific, we shall consider two gluons with momenta qi and (72, which are rescaled by a 
factor A 

qi -> Xqi , 92 M2 , (B.l) 



and we will study the limit A — > 0. As explained in [57|, the matrix element factorizes as follows 



(ai, 02; /ii,/i2|A^g,g,ci,...,c„(gi, (72,^1, •••,Pn)) 9'^p'^''J^l';ll{qi,q2)\Mcu...,cApi, ■■■,Pn)) , (B.2) 

where g is the strong coupling constant, fi the dimension unit in dimensional regularization (intro- 
duced through the explicit dependence in the coupling constant), and the two-gluon soft current 
J^lZ^1i^<l2) is given by 

j^Z'(quq2) = (.r), (.2)} + ./.....3 E rr I J;^^^ . ^^^^ 



Pt ■ (qi - 92) 



{pi ■ qi){pi ■ q2) qi-q2 



(B.3) 



2 [pi ■ {qi + q2) 

with the eikonal current defined as 

J^(g)=ET,:^. (B.4) 

The algebra of the colour operators has been discussed at length in ^ . The expression Eq. (jB.3p 
is as in jST] and can be derived by taking into account all diagrams with the two soft gluons attached 
to a hard parton line through an eikonal coupling. The triple gluon vertex has to be treated exactly, 
since all momenta are of the same order. The chosen class of diagrams is shown to be sufficient 
by power counting in a physical gauge. Moreover, contraction with physical polarization vectors 
has been used to eliminate terms proportional to the soft gluon momentum. 

The difference between massive and massless cases occurs, when squaring the matrix element. 
The factorization formula then contains the factor 

[j;r(gi,92)]^d^"(gi)d''''(g2)J."r(9i'92) = ^ {J'(gi), J'(g2)} -C^ E ^.-T, 5,, (91,92) + ... , 

(B.5) 

where d'^'^{q) is the polarization tensor obtained by summing over gluon polarizations. Due to 
current conservation, we can make the replacement d'^'^{q) — >■ — g'"^. The terms vanishing when 
acting on a physical matrix element are denoted by the dots at the end of the above equation. 

In order to recast what is essentially the square of the two-gluon current in Eq. (jB.3p into the 
form of the right hand side of Eq. (IB. 51) . some colour algebra is needed (rightfully called "quite 
cumbersome" by the authors of 57]). The process is simplified substantially by the use of the 
following two identities 

in"'"' [{T^^\T^'},T;!'] - 2CAT,-T,{d,k-Sjk); (B.6) 
{{Tt\T;'-}An\Tr}} + {{Tr,Tr},{T^\T;^-}} 

= 8 {Ti ■Tk,Tj • T;} + 2Ca SSuSjkTi ■ Tj + 35^4; • 
—2Sij5jkTi ■ Ti — 26ijSjiTi ■ Tk — 2{SikSki + SjkSki)Ti ■ Tj 
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The result for the Sij{qi,q2) function can be spUt into two parts 



Sij{qi,q2) ^S"l °((?i,(?2) + (m'^ (72) + ^"^"(qi, 92)) 



(B.7) 



where the first term has already been given in [57| and reads 

(1 - e) Pi - qi pj ■ q2+Pi- q2 Pj ■ qi 



5™=° (91 '92) 



[qi ■ 92)^ p^ ■ {qi + 92) Pj ■ [qi + 92) 



(Pi-Pj) 



2pi- qi Pj ■ q2Pi- <72 Pj ■ qi 



Pi ■ qi Pj • 92 + Pi • 92 Pj ■ qi 



Pi ■ Pj 
2 91 ■ 92 



Pz ■ (91 + 92) Pj ■ (91 + 92) 
2 1 



Pi ■ 91 Pj ■ 92 Pj -qiPi- 92 Pi ■ (91 + 92) Pj ' (9i + 92) 



X 4 + 



{Pi ■ 91 Pj ■ 92 +Pi ■ 92 Pj ■ qif 

Pi ■ 91 Pj ■q2Pf 92 Pj ■ 91 



(B.6 



The second contribution in Eq. (jB.7l) is new and represents the additional terms generated by 
non-vanishing masses. The relevant function is 



5l"'^°(9i,92) = 



Pi ■ Pj Pj ■ (91 + 92) 



4 91 ■ 92 Pi • 91 ■ 92 '2. Pi- 9i Pj • 92 Pi ■ 92 Pj ■ 9i Pi ' (9i + 92) 



2 91 • 92 Pi • (91 + 92) Pj ■ (91 + 92) VPi • 91 Pj -92 Pi- 92 Pj • 9i 



{Pj -qi? 



(Pj -92)^ 



(B.9) 



Appendix C. Collinear behavior in the double-soft Umit 

Due to the particular phase space decomposition introduced in the singular matrix element 
limits that need to be considered in the construction of the subtraction terms, are covered directly 
by the formulae from [Ht} (aside from the modification given in the previous appendix for the case 
of massive partons). Nevertheless, one case turns out to be slightly inconvenient. Indeed, the 
double-soft limit followed by the collinear limit of the two soft partons, although obtainable with 
the formulae of |Appendix B[ requires a careful evaluation, because of the presence of an apparent 
quadratic divergence ~ l/(9i • 92)^ in Eq. ()B.8| . Of course, the actual leading divergence is only 
logarithmic as can be checked using colour conservation. To avoid unnecessary complications, we 
propose to use an iterated limit in which the partons become collinear first, and then produce a 
soft gluon, which interacts via the usual eikonal current. This is justified by colour coherence of 
soft emission in the collinear limit, which is usually exploited to derive the soft-collinear limit in 
which a pair of partons become collinear, and a gluon, not belonging to the pair, becomes soft. 
The result for our case can be written as follows 

2 

|A^ci,C2,ai,...,a„(9l,92,Pl, ■•■,P«)P - 9 V^' (-^ai , . . . ,a„ (pi , ■ . ■ , P«) | (C.l) 

Sl2 

(Pl, ■•■,Pn)) , 

where C1C2 — qq or gg, S12 = (91 + 92)^, is the eikonal current defined in Eq. (IB. 41) . and 
Pj!"^^ {z, k±; e) is the d-dimensional polarized Altarelli-Parisi splitting function given in Eqs. (|A.6[ 

roitaroi). 

The factorization formula demonstrates the usual spin correlations, which are transferred here 
to the eikonal currents and not directly to the matrix element, since the nearly on-shell gluon 
is fully described by the current in the soft limit. One might wonder why the spin correlations 
survive the soft limit, since they do not at the next-to-leading order. The reason is that the 



40 



double-soft limit cannot be defined with the momenta of the collinear pair alone, because the 
splitting functions depend only on the ratio of the energies of the two partons. Therefore, as long 
as this ratio remains constant, from the point of view of the collinear limit, we are considering two 
hard partons. 

Appendix D. Born level cross sections for top quark pair production 

Although text book material, we reproduce these cross sections here for convenience of the 
reader. We have 



with 



IT TpCp 
6 N 

TT Tp 



12 

Ca 



Pp (2 + p) , 



(4 + 4p-2p^) - ln^-4-4p 



P 1-/3 



(D.l) 



(D.2) 



(D.3) 



and p = 1 - . 

Appendix E. Software 

The results obtained for the present publication have required the use of numerous software 
systems. We list them here 

• DiaGen /IdSolver, our own private system for diagram generation, analysis and evaluation, 
has been used for the generation of the cut diagrams and reduction of the integrals needed 
to compute the volume of the phase space; 

• Fermat ^6§j , an algebra system, is the rational function algebra library of DiaGen /IdSolver, 
and has been used in the reduction of the phase space integrals; 

• Form [g^, has been used for the algebraic simplification of the diagrams, mostly Dirac 
algebra and color factor evaluation, for which the package Color. H has proven useful; 

• FormCalc [t^, the backbone of FeynArts [7i|, has been used for the low level formatting 
of Fortran code generated by Mathematica; 



• Helac/Phegas [61|, |62|, [72], has been used for tests of the matrix elements at specified 
phase space points and checks of the numerical integration routines; 

• Mathematica, has been used for the derivation of the subtraction and integrated subtrac- 
tion terms, convergence tests with very high numerical precision, and generation of Fortran 
code; 

• Intel Fortran Compiler, although not essential for this project, we have used its quadru- 
ple precision functionality to spare some minor effort in implementing interfaces to external 
libraries (see comments in Section |4|); 

• Parni (g^ I adaptive Monte Carlo random number generation optimizer, has been used in 
the numerical integration routines; 



• Ranlux [73|,lIJ], the classic random number generator. 
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